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Abstract 


This  paper  presents  two  algorithms  for  solving  sparse  nonlinear  systems  of  equations:  the 
CM-successive  column  correction  algorithm  and  the  modified  CM-successive  column  correction 
algorithm.  A  g-superlinear  convergence  theorem  and  an  r-convergence  order  estimate  are  given 
for  both  algorithms.  The  numerical  results  indicate  that  these  two  algorithms,  especially  the 
modified  algorithm  are  probably  more  efficient  than  some  currently  used  algorithms. 


1.  Introduction. 


Consider  a  nonlinear  system  of  equations 

F(x)  =  0,  (1-1) 

where  F:  R*-*R*  is  continuously  differentiable  on  an  open  convex  set  DC.R* ,  and  the  Jacobian 

matrix  F  '(x)  is  sparse.  To  solve  the  system,  the  following  iteration  is  considered: 

xk+i  =  xk  _  Bt'Fiz*)  ,  *=0,1,...,  (1.2) 

where  B*  is  an  approximation  to  the  Jacobian  with  the  same  sparsity  structure. 

For  convenience,  we  rewrite  (1.2)  as 

x  =  x-B~lF(x)  ,  (1.3) 

where  x  and  x  indicate  the  current  iterate  and  the  new  iterate  respectively,  and  B  is  an  approxi¬ 
mation  to  the  Jacobian. 

Currently,  there  are  several  algorithms  to  get  a  sparse  approximation  to  the  Jacobian.  In 
this  paper  we  will  discuss  three  types  of  algorithms. 

(1)  Schubert’s  algorithm.  In  1970  Schubert  [17]  gave  a  sparse  modification  of  Broyden’s 
update.  Broyden  [2]  also  gave  this  algorithm  independently.  In  order  to  present  Schubert’s  algo¬ 
rithm,  we  introduce  the  following  notation  concerning  the  sparsity  pattern  of  the  Jacobian: 

Definition  1.1.  For  j= 1,2,...,  define  the  subspace  Z,C.Rn  determined  by  the  sparse  pattern  of  the 
j th  row  of  the  Jacobian: 

Zj^{vERu:  e,rv=0  for  all  i  such  that  [F'(x)]y,=0  for  all  xER*}, 
where  e,-  is  the  tth  column  of  the  n  X  n  identity  matrix.  Define  the  set  of  matrices  Z  that  preserve 

the  sparsity  pattern  of  the  Jacobian: 

Z  =  {AEL(Ru):  At CjEZj  for  j  =  l,2,...,n  }. 

Definition  1.2.  For  j= l,2,...,n,  define  the  projection  operator,  DjEL{R *),  that  maps  i?"  onto  Zf. 

Dj  —  diag  •  •  •  i  d/»), 

where 


1,  if  e,-  E  Zjt 
0,  otherwise. 


For  a  scalar  a ER ,  define  the  pseudo-inverse: 


a  = 


if  a  7^  0, 
if  a  =  0. 


Now  Schubert’s  update  can  be  written  as 


B  =  B  +  £ ( [* ] /l*  1  y )+  c/( »  -  Bo )[«]/, 

i-i 

where  [s]y  —  Z>y«,  «  =  z  -  x  and  y  =  F(x)  -  F(x). 


(1-4) 


Let 


<3.it  =  {A6L(f?*):  >lu=»,  for  vectors  u,  vER*}. 

The  following  theorem,  which  we  will  use  later,  was  proved  by  Reid  [15]  and  Marwil  [9]  indepen¬ 
dently. 


Theorem  1.1.  Let  BEZ ;  y ,  sERk  with  sj£ 0.  Define  B  by  (1.4).  Then  B  is  the  unique  solution  to 

wa{\\6-B\\,:6eQ,,.nZ}t  (1.5) 

where  1 1 . 1 1  p  indicates  the  Frobenius  norm  of  a  matrix. 

The  advantage  of  Schubert’s  algorithm  is  that  at  each  iteration  only  one  function  value  is 
required  and  it  is  g-superlinearly  convergent  (see  Marwil  [9]).  However,  it  frequently  requires  more 
iterations  than  finite  difference  algorithms.  Moreover,  it  may  not  be  a  good  approximation  to  the 
Jacobian  when  the  problem  is  badly  nonlinear,  especially  when  the  current  step  is  far  away  from 
the  solution.  Therefore,  p*  =  -B^1F(zi)  may  not  be  a  descent  direction  of  the  functional 

/(*)»— 1|  F(*)  || 2,  where  ||.||  denotes  the  /2  vector  norm.  In  this  case,  it  may  be  not  good  to 
2 

use  a  line  search  with  Schubert’s  algorithm. 

(2).  Finite  difference  algorithms.  In  general,  a  finite  difference  algorithm  can  be  formulated 
as  follows:  obtain  direction  vectors  dl,d2,  .  .  .  ,df  such  that  B  can  be  determined  uniquely  by  the 
equations 

Bdi  =  F{x+di)  -  F{x),  i=l,2,...,p. 
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In  this  paper,  we  assume  that  it  is  not  convenient  to  evaluate  the  function  values  element 
by  element,  instead  we  only  evaluate  the  value  of  F(x)  as  a  single  entity.  This  is  reasonable  since 
in  practice  it  is  very  common  that  the  components  of  F(x)  have  expensive  common  sub¬ 
expressions.  In  this  case,  to  reduce  the  number  of  function  evaluations,  Curtis,  Powell,  and  Reid 
[4]  proposed  a  finite  difference  algorithm,  called  the  CPR  algorithm,  which  is  based  on  a  partition 
of  the  columns  of  the  Jacobian.  Coleman  and  More  [3]  associate  the  partition  problem  with  a 
graph  coloring  problem  and  gave  some  partitioning  algorithms  which  can  make  the  number  of  the 
function  evaluations  optimal  or  nearly  optimal. 

Following  Coleman  and  More',  we  give  some  definitions  concerning  a  partition  of  the 
columns  of  the  Jacobian. 


Definition  1.8.  A  partition  of  the  columns  of  a  matrix  B  is  a  division  of  the  columns  into  groups 
Ci,c2,...,c,  such  that  each  column  belongs  to  one  and  only  one  group. 


Definition  l.j.  A  partition  of  the  columns  of  a  matrix  B  is  consistent  with  the  direct  determina¬ 
tion  of  B  if  whenever  6,-y  is  a  nonzero  element  of  B,  then  the  group  containing  column  j  has  no 
other  column  with  a  nonzero  element  in  row  * . 

As  an  example  we  consider  the  tridiagonal  structure 
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(1.6) 


A  consistent  partition  of  the  columns  of  the  matrix  is  Cj  =  {1,  4},  c2  =  {2,  5},  and  c3  =  {3,  6}. 

The  CPR  algorithm  now  can  be  formulated  as  follows:  for  a  given  consistent  partition  of  the 
columns  of  the  Jacobian,  obtain  vectors  dud2,.-,dr  such  that  B  is  determined  uniquely  by  the 
equations 
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Bdt  =  F(x+d{)  -  F(x)  a  .=1,2,. ..,p  .  (1.7) 

Notice  that  for  the  CPR  algorithm,  the  number  of  function  evaluations  at  each  iteration  is  p  +  1. 

Since  the  partition  of  the  columns  of  the  Jacobian  plays  an  important  role  in  the  CPR  algorithm, 

we  call  the  CPR  algorithm  based  on  Coleman  and  Morels  algorithms  the  CPR-CM  algorithm. 

For  the  consistent  partition  given  in  example  (1.6),  if  we  take 

dt=(A  ,0,0,  A  ,0,0)r, 
d2  =  (  0  ,  h  ,  0 , 0  ,  h  ,  0  )T, 
d3  =  (  0 , 0  ,  h  ,  0 , 0  ,  h  )T, 

then  B  is  determined  uniquely  and  the  number  of  function  evaluations  required  at  each  iteration 
is  4. 

The  advantage  of  the  CPR  algorithm  is  that  it  usually  requires  fewer  iterations  than 
Schubert’s  algorithm.  However,  it  requires  more  function  values  at  each  iteration  than  Schubert’s 
algorithm. 

(3).  The  successive  column  correction  algorithms. 

Polak  [13]  gave  a  successive  column  correction  algorithm  for  unconstrained  minimization. 
Feng  and  Li  [7]  developed  a  successive  column  correction  algorithm  for  nonlinear  system  of  equar 
tions,  which  is  called  the  column-update  quasi-Newton  method.  Using  this  algorithm,  columns  of 
Bk  are  displaced  by  differences  successively  and  periodically.  At  each  iteration,  only  two  function 
values  are  required,  but  only  one  column  is  displaced. 

In  this  paper,  we  propose  two  algorithms:  the  CM-successive  column  correction  algorithm 
and  the  modified  CM-successive  column  correction  algorithm.  The  former  is  based  on  Coleman 
and  Morn’s  algorithm  and  the  column-update  algorithm.  The  latter  is  a  combination  of  the  CM- 
successive  column  correction  algorithm  and  Schubert’s  algorithm.  Both  algorithms  require  only 
two  function  values  at  each  iterative  step.  Our  numerical  results  show  that  the  CM-successive 
column  correction  algorithms,  especially  the  modified  one,  are  probably  more  efficient  than  the 
CPR  algorithm  and  Schubert’s  algorithm. 
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The  CM-successive  column  correction  algorithm  is  given  in  Section  2.  A  Kantorovich-type 
analysis  for  this  algorithm  is  given  in  Section  3.  A  g-superlinear  convergence  result  and  an  r- 
convergence  order  estimate  of  the  CM-successive  column  correction  algorithm  are  given  in  Section 
4.  The  modified  CM-successive  column  correction  algorithm  is  given  in  Section  5.  Some  numeri¬ 
cal  results  are  given  in  Section  6. 

In  this  paper,  for  a  sparse  matrix  B,  we  use  M  to  denote  the  set  of  pairs  of  indices  («,  j), 
where  b,-}-  is  a  structurally  nonzero  element  of  B,  i.e. 

M  hi*  0}. 


2.  The  CM-Successive  Column  Correction  Algorithm  and  its  Properties. 

Given  a  consistent  partition  of  the  columns  of  the  Jacobian,  which  divides  the  set 
{1,  2,  ...,  n}  into  p  subsets  ct,  e2,  epi  let 


where 


(2.1) 


«*=*  (  mod  p  ),  *==1,2,..., 

and  let 

yk  =  F{xk  -(-  dk)  -  F(xk).  (2.2) 

The  CM-successive  column  correction  algorithm  can  be  formulated  as  follows:  If  k<p,  then  for 

jecit  the  jth  column  of  Bk  is  determined  uniquely  by  the  equation 

Bkdk  =  yk,  (2.3) 

and  the  other  columns  of  B*  are  equal  to  the  corresponding  columns  of  B*_ j.  If  k  >p ,  the  columns 
of  Bi,  are  displaced  as  described  above  successively  and  periodically.  In  other  words,  for  the 

jth  column  of  B^  is  determined  uniquely  by  (2.3),  and  the  other  columns  of  B*  are  equal  to  the 
corresponding  columns  of  B *_t. 

For  example  (1.6),  at  the  first  iteration  we  displace  the  first  group  Cj  =  {1,  4}.  At  the 
second  iteration  we  displace  the  second  group  e2  =  (2,  4}.  At  the  third  iteration  we  displace  the 
third  group  Cj  =  {3,  6},  and  then  we  displace  the  three  groups  successively  and  periodically. 
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The  CM-successive  column  correction  algorithm  with  a  global  strategy  is  given  below. 


Algorithm  2.1.  Given  a  consistent  partition  of  the  columns  of  the  Jacobian,  which  divides  the  set 
{l,  2,  ...,  n}  into  p  subsets  Cj,  e2,  — ,  cf  (for  convenience,  c,-,  »=l,2,...,p,  indicates  both  the  sets 
of  the  columns  and  the  sets  of  the  indices  of  these  columns),  and  given  an  x°€R *  and  a  nonsingu¬ 
lar  matrix  B0,  which  has  the  same  sparsity  as  the  Jacobian,  at  the  initial  step: 

(1) .  Set  /  =  0. 

(2) .  Solve  B0s°  =  -F(x°). 

(3) .  Choose  xl  by  x1  =  x°  -I-  s°,  or  by  a  global  strategy. 

At  each  iteration  it  >  0: 

(1) .  Choose  a  scalar  A* . 

(2) .  If  /  <  p,  then  set  1  =  1  +  1,  otherwise  set  /  =  1. 

(3) .  Set 

dk  =  £A*ey. 

(4) .  If  j£ct  and  (i,  /)  6  M,  then  set 

b £  =  rW(F(x*  +  d")  -  F(xk)),  (2.4) 

rl 

otherwise  set 

A* .  —  A*.-i 

vtj  utj  > 

where  Bk  =  [6*]. 

(5) .  Solve  —  -F{xh). 

(0).  Choose  xk+l  by  x*+1  =  x*  +  «*,  or  by  a  global  strategy. 

(7).  Check  for  convergence. 


Let 


Then 


J>  =  f  F'(x"+tdk)dt  . 


(2.5) 


0 


Jkdk  =  yk  .  (2.6) 

Let  Jk  =  -  Since  Jk  has  the  same  sparsity  as  the  Jacobian,  by  (2.6),  we  have  that  if 


(/,  m)  €  M,  then 


JL  = 


T  k 

V 


(2.7) 


where  m6c,(.  Comparing  (2.7)  with  (2.4),  we  have 


Bk  tj  —  Jk  Cj  , 

for  jGe.y 

The  CM-successive  column  correction  algorithm  is  also  an  update  algorithm,  and  the  update 
can  be  written  as: 


je% 

From  (2.8),  it  is  easy  to  get  the  following  result: 


Bk  =  BUI  -  2  «/*/)  +  E 


(2.8) 


Lemma  2.2.  Let  Bk,  Jfc=l,2,...,  be  generated  by  Algorithm  2.1.  If  k>p ,  then 


Bk  —  E  E  h  et  eiT 

/=*-?+ 1  let, 


(2.9) 


To  study  the  properties  of  our  algorithms,  sometimes  we  assume  that  F'  satisfies  the  fol¬ 
lowing  Lipschitz  condition:  there  exist  a,->0,  i'=l,2,...,n  such  that 


||(F'(x)  -  F'(y))e,- 1|  <a{  ||x  -  y  ||,  x,yED. 


Let  a==(E“?)2.  Then,  it  follows  from  (2.10)  that 
1  —  1 


|F'(x)  -  F'(y)||f  <a||*  -  y  ||,  x,yED. 


(2.10) 


(2.11) 


Theorem  2.3.  Let  F'  satisfy  Lipschitz  condition  (2.10).  Also  let  {xy)y=0C D  and  let  {By}*=0  be 

generated  by  Algorithm  2.1  with  |  A*  |  <  — ^  1 1  x*  —  x*-1 1 1 .  If 

v  « 

{xJ+«f,}*_iC.D,  then  for  k>p, 


II  Bk-F'(xk)  ||  f  <a  E  ll*y-*y_1| 

/— k-p+i 


(2.12) 
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Proof.  By  (2.5),  (2.1)  and  Lipschitz  condition  (2.10), 


||(F'(x») -/„)*,•  || 

=  ||(/0V'(afm+ldm)-F'(zm))dO«yll 

<aJ0,|Mm||^  =  -f-|irf"H  (213) 

=  ^-11  E  *m«yll 

<  ?L^\hm\  <  atj  \  \xm  -  xm~l  1 1  , 

where  k-p+ 1  <  m  <  Jfc.  It  follows  from  (2.9)  and  (2.13)  that 
\\F'(x>)-Bt\\t 

=  E  ||  E(^V)-B*Ke/||| 

m=*-p+i  ye«, 

=  E  E  ||(F'(x*)-/m)ey||2 

m=*-p  +  l  ject" 

<  E  E  (ll(*,'(**)-*,'(*m))«yll  +  H(F'(x"Km)ey||)2 

m— *-p+l  /€«,■ 

<  E  E  «J(II*‘  -*"  ||  +  ||*m  -xm-l||)2  (214) 

m=*-p+i  ye<!lm 

<  E  E  «y(  E  ll*'-*'-'ll)2 

m— *-p+l  ;€«,  J=*-p+l 

=  a2(  £  ll*‘  -xM||)2- 

i>.t-p+i 

Then,  (2.12)  follows  from  (2.14). 

To  start  iteration  (1.2)  for  a  given  x°£D ,  an  initial  matrix  B0  is  needed.  We  suggest  using 
the  CPR-CM  algorithm  to  get  B0  since  it  is  easy  to  implement  after  we  have  a  consistent  parti¬ 
tion  of  the  columns  of  the  Jacobian. 
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3.  A  Kantorovich-Type  Analysis. 

By  means  of  Theorem  2.3,  we  have  the  following  Kantorovich-type  analysis  for  the  CM- 
successive  column  correction  algorithm. 


Theorem  8.1.  Assume  that  F'(x)  satisfies  Lipschitz  condition  (2.10).  Let  x°£D,  and  let  B0  be  a 
nonsingular  n  X  n  matrix  such  that 


and 


II 50  -  F  '(z°)  ||  f  <  S,  \\B?\\F<0,  ||5o1F(x°)||  <  f?  , 

l  =  a  Pv  <  2_ 

(1  -  Z08f  “  6  ’ 


(3.1) 


If  S(x°,  21*)  C  D ,  where 


then  {xk },  generated  by  the  CM-successive  column  correction  algorithm  with 

2  # 

|  A*  |  <— ||  **  -  x*_1 1|  and  without  any  global  strategy,  converges  to  i*,  which  is  the  unique 

vn 

root  of  F(x)  in  5(x°, 7)r\D,  where 


t  = 


1-06 

a  P 


1  + 


2  a0rf 


Proof.  Consider  the  scalar  iteration 

tk+i-tk  t0  =  0,  k  —0,1,2,  ••  •  ,  (3.3) 

where 

/( =  +  |  •  (3.4) 

It  is  easy  to  show  that  the  sequence  {<*  }  satisfies  the  difference  equation 

**+i-**  =  $0  [—  (^*  —  **-i) « <*-i -b  (**-*»- i)  »  *=1,  2,  •  (3.5) 

From  this  equation,  it  follows  that  {t*}  is  a  monotonic  ally  increasing  sequence  and 
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lim  =  t*  , 

k  -*oo 


where  t '  is  the  smallest  root  of  / (/). 

Now,  by  induction,  we  will  prove  that 

||z*+1  —  x*  ||  <  ,  4  =1,2,  •  • 

{**}  C  S(x°,t‘),  k  =  1,2,  ■  ■  ■  , 

{xk  +  dk}  C  S(x°, 2t‘), 

and 

||5*~1|l  <  3/9,  4  =  1,2, 

For  1  =0,  we  have 


(3.6) 

(3.7) 

(3.8) 

(3.9) 


II  <  n  =  ti-fo  <  • 


Thus, 


l|x1  +  d1-x°||  <  ||«1-«°||  +  |l'l|l  <2||x1-x°||  <  2f*  . 
Suppose  (3.6)  holds  for  k  =0, 1, ...,  m  -1.  Then, 


ll*"-*°ll  <  E  <  f. 

1=0 

Therefore,  xm  £S  (x°,  <*),  and 

{xm  +  dm}  C  5(x°,2t')  . 

From  the  proof  of  Theorem  2.3,  it  can  be  seen  that  for  all  k, 


l|S*-F'(x*)||/.<||5a-F'(x°)||f  +  «£  llx'-x'-Ml 

i'=  o 


Therefore, 


||V(^-B0)|| 

<  1 1  So'1  ||  f(  II  Sm -F '(x")  ||  f  +  ||F'(x’")-F'(x0)||f  +  |  |F '(x°)-F0 1|  f ) 

<  #2a  1 1  x,+1  -  x*  1 1  +  26) 

i-0 

<0{2ar  +  26)<^)  =  l. 


(3.10) 
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Thus,  by  Dennis  and  Schnabel’s  Theorem  3.1.4  [6,  p.45], 


l|B"  11  -  T^/3  ~3^  ' 

Hence, 

||i"+1-*"  || 

<  ||  B*1  ||f  ||F(*")-F(*m-1)-5m_i(*m  || 

<  3(3[^||x"-x"-1||+aE  1 1  xi+i-x'  ||  +  6]  ||  arm  -ar'"-1 1 1 

A  »=o 

—  3/?[— (tm  -  lm_i)  +  a  lm_i  +  6](fm  -  *m-l)  =  Ln-t-l- tjn  ■ 

This  completes  the  induction  step.  By  (3.6),  it  is  easy  to  show  that  there  is  an  x*  E D  such  that 

lim  x*  =  x  *  . 

k-*oo 

The  uniqueness  of  x*  in  S(x°,t  )f\D  can  be  obtained  from  Ortega  and  Rheinboldt’s  Theorem 
12.6.4  [12,  p.425]  by  setting  A(x)  =  B0. 

4.  Local  Convergence  Properties. 

To  study  the  local  convergence  of  our  algorithms,  we  assume  that  F.D  C  R*  —+R*  has  the 
following  property: 

There  is  an  x‘  ED  ,  such  that  F(x*)  =  0  and  F ' (x*)  is  nonsingular.  (4.1) 

Theorem  j.l.  Let  F:D  C  R*  -+ Rn  satisfy  (4.1),  and  let  F"  satisfy  Lipschits  condition  (2.10). 

2 

Also  let  {x*}  be  generated  by  Algorithm  2.1  with  |  A*  |  <— — -||x*  -  x*_l  1 1  and  without  any 

vn 

global  strategy:  Then,  there  exist  «,  S  >  0  such  that  if  x°ED  and  B0  satisfy 

II  z°_x*  ||  <  « ,  ||B0-F'(x*)||f  <6, 

then  {x*}  is  well  defined  and  converges  g-superlinearly  to  x*. 

Proof.  Notice  that  when  e  and  6  are  small  enough,  we  have  that  A  <  — ,  f}6  <  —  and  that 

6  3 

S(x°,2t*)  C  D  where  A,  0  and  1*  are  defined  in  Theorem  3.1.  Therefore,  by  Theorem  3.1, 
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x*  +  dk  E  D  ,  k=0,  1,  •  •  •  . 

By  (2.8), 

Bk-F'(a t*) 

•(,■))(/-  E «,«,’)  +  Et4-F'(*'IK‘'  <4'2> 

ye<,t  j€«,t 

Thus, 


<a(||x*-x‘||  +±||d*|l)  (4-3) 

<«(ll**-*'ll  +  ll**-*‘-1H) 

<a(2||*‘-«*||  +  ||*‘-1-x*||). 

Let  <r(x*_1  ,  x*)  =  max  {  1 1  x*  -  x  *  1 1  ,  1 1  x*-1  -  x‘  \  |  }.  Then  it  follows  from  (4.2)  and  (4.3)  that 

||fl»-F'(x*)llr  <  ||B*-1-F'(x')||,.+  ||/*-F'(*')llr 
<  -  F'(x*)||f  +  3aa(x*_1  ,  x*)  . 

Thus,  by  Dennis  and  Morels  [5]  Theorem  5.1,  we  know  that  {x*}  converges  at  least  ^-linearly  to 
* 

X  . 

According  to  Dennis  and  Morels  [5]  Theorem  3.1,  to  get  <7 -super]  inear  convergence,  we  need 
only  to  prove  that 


From  (2.12),  it  follows  that 


lim 


11  (#*  -■F'(x*))(x*+1-x*) 


=  0  . 


lim  ||B*  -F'(x*)||=0. 

•  — >oo 

This  implies  (4.4). 


(4.4) 


(4.5) 


Theorem  4-8.  Assume  that  F  satisfies  the  hypotheses  in  Theorem  4.1.  Then  the  r-convergence 
order  of  Algorithm  2.1  is  not  less  than  r,  where  r  is  the  unique  positive  root  of 

=  o  . 


la 


Proof.  Notice  that  (4.5)  implies  that  there  exist  k0  and  f3  >  0  such  that  1 1  B*  1 1 1  <  ft  for  all 
k  >  k0.  Thus,  by  Theorem  2.3, 

||  x*+1-ar*  ||  =  ||**-z*-B*-1E(**)ll 

+  (l|F'(0-FV)ll^+II^V)-ftllrfll**-*'ll} 

<  a  ll**-*'ll  +  «  E  ||*/+1-*y||}||**-*'ll 

1  /=*-* 

<  |«/?(  E  ll*y-*'ll)ll**-**ll  • 

1  J=k-p 

Thus,  the  desired  result  follows  from  Ortega  and  Rheinboldt’s  Theorem  9.2.9  [12,  p.29l]. 

5.  The  Modified  CM-Successive  Column  Correction  Algorithm. 

Estimate  (2.12)  shows  that  when  p  is  small,  5*  is  a  good  approximation  to  F'(i :*).  How¬ 
ever,  B/,  still  retains  information  from  the  previous  p  steps.  Therefore,  the  following  question  is 
reasonable:  Can  we  have  a  better  approximation  to  F  '(z*)  without  more  function  evaluations? 
Notice  that  when  we  get  Bk  by  Algorithm  2.1,  we  did  not  use  the  value  of  F(ii).  The  main  idea 
of  the  modified  CM-successive  column  correction  algorithm  stated  below  is  to  use  all  the  informar 
tion  we  already  have  to  improve  our  approximation  to  F'{xk). 

Algorithm  5.1.  Given  a  consistent  partition  of  the  columns  of  the  Jacobian,  a  vector  x°  and  a  non¬ 
singular  matrix  B0  with  the  same  sparsity  as  the  Jacobian,  at  the  initial  step: 

(1) .  Set  1=0  and  B0  =  B0. 

(2) .  Solve  B0s°  =  -F(x°). 

(3) .  Choose  x1  by  x1  =  x°+»°,  or  by  a  global  strategy. 

At  each  iteration  Jfc  >  0: 

(1) .  Update  by  Algorithm  2.1  to  get  B*. 

(2) .  Update  B*  by  Schubert’s  update  to  get  B*. 

(3) .  Solve  B*«*  =  -F(xk). 
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(4) .  Choose  x*+1  by  x*+1  =  x*+a*,  or  by  a  global  strategy. 

(5) .  Check  for  convergence. 

Our  numerical  results  show  that  Algorithm  5.1  usually  requires  fewer  iterations  than  Algo¬ 
rithm  2.1.  Especially,  when  the  problem  is  not  well  behaved,  and  a  global  strategy  is  used,  the 
modified  algorithm  behaves  significantly  better  than  Algorithm  2.1.  The  cost  of  the  improvement 
is  the  computation  of  Schubert’s  update.  However,  since  the  Jacobian  is  sparse,  Schubert’s 
update  requires  only  0(n)  operations.  We  feel  that  it  is  worth  doing  this  rather  than  computing 
more  function  values  and  solving  more  linear  systems. 

Now  we  will  briefly  discuss  the  convergence  properties  of  Algorithm  5.1.  Let 

(5.1) 

Since  /*  performs  exactly  the  same  as  the  secant  factor  — ‘  J  \  L  in  one  dimensional  prob- 

x  -x*  1 

lems,  we  call  Jk  the  secant  operator.  It  is  easy  to  show  the  following  result. 

Lemma  5.1.  If  F'  satisfies  Lipschitz  condition  (2.10),  then 

l|J.-f'(«‘)llF<f  •  (5-2) 

Estimate  (5.2)  shows  that  /*  is  a  good  approximation  to  F  '(x*)  when  1 1  x*  —  xil  ||  is  small. 

Theorem  5.2.  Let  F'  satisfy  Lipschitz  condition  (2.10).  If  {5*}  and  {Bk  }  are  generated  by  Algo¬ 
rithm  5.1,  then 

11^-7*11/-  <  ||Bt-Tt|U  (5-3) 

If,  in  addition,  fl*  B/,,  then  the  strict  inequality  holds. 

Proof.  Since  7, i€Qt,,C\Z,  by  Theorem  1.1,  we  have 

1 1  Bk  -  *4  1 1  *  +  1 1  Bk  -  Bh  1 1  p  =  1 1  Bk  -  /*  1 1  p .  (5.4) 

Then,  (5.3)  follows  from  (5.4). 

Notice  that  in  general,  Bk  Bk.  Therefore,  by  Theorem  5.2,  5*  is  usually  closer  to  the 
secant  operator  J *  than  Bk.  Thus,  fl*  should  be  a  better  approximation  to  the  Jacobian  than  5* 
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when  Bit  retains  some  information  from  previous  steps.  But  theoretically,  we  can  not  get  a  better 
estimate  for  ||B*  -  F  '(x*)  1 1  f  than  that  for  ||B*  -  F  '(x*)  ||  However,  we  can  get  the  follow¬ 
ing  result: 

Theorem  5.8.  Let  F  :  R *  — ►  R *  satisfy  Lipschitz  condition  (2.10).  Also  let  {B*}  and  {a:*}  be 
generated  by  Algorithm  5.1.  Then, 

\\Bk-F'(xk)\\F<2a  £  ||*y-*/-1||  •  (5.5) 

j=k-p+i 

Proof.  By  (5.3), 

\\Bk-F'(xk)\\F 

<  ||*- J*||,  +  ||7t-fV)ll, 

<  ||fl*-7*||,+  ||7fc-FV)ll/- 

<  ||5*-fV)II,  +  2||7*-fV)II,. 

Then,  from  (2.12)  and  (5.2),  we  obtain  (5.5). 

From  estimate  (5.5),  it  is  easy  to  prove  that  Algorithm  5.1  has  at  least  the  same  local  con¬ 
vergence  properties  as  Algorithm  2.1. 

8.  Numerical  Results. 

We  computed  some  examples  with  tridiagonal  Jacobians  by  the  CPR  algorithm,  Schubert’s 
algorithm,  Algorithm  2.1,  and  Algorithm  5.1.  In  this  section,  we  compare  the  numerical  results 
from  these  four  algorithms.  The  global  strategy  we  used  in  computing  the  examples  is  the  line 
search  with  backtracking  strategy  (see  Dennis  and  Schnabel  [8,  p.128]).  For  the  CPR  algorithm,  if 
p*  =  -BTlF(i t*)  is  not  a  descent  direction,  then  we  try  -p*.  If  it  is  not  a  descent  direction  either, 
then  the  algorithm  fails.  For  the  other  algorithms,  if  p*  is  not  a  descent  direction,  then  we  try 
-p*.  If  it  is  not  a  descent  direction  either,  then  we  try  the  CPR  direction.  If  the  CPR  direction 
fails,  then  the  algorithm  fails.  In  the  CPR  algorithm,  Algorithm  2.1  and  Algorithm  5.1,  at  step  k, 
we  use  different  A*  for  each  component  of  x*  instead  of  one  uniform  A*.  According  to  Dennis  and 
Schnabel  [6,  p.98],  we  choose 
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The  stopping  test  we  used  is 


h *  =  V machepa  x* . 


,k+l 


max 


~rfl 


!<*'<»  max{  |  x*+1  | ,  typxi} 
and  we  choose  c  =  1CT6.  We  used  double  precision,  and  the  machine  precision  is  2.22</-18. 


Example  6.1  was  given  by  Guangye  Li  [8],  and  it  can  be  seen  to  be  an  extension  of  the 
Rosenbrock  [16]  function  (also  see  More',  Garbow  and  Hillstrom  [ll])  to  nonlinear  system  of  equar 
tions  with  tridiagonal  structure.  Example  6.2  was  given  by  Broyden  [l]  (also  see  More',  Garbow 
and  Hillstrom  [11]).  Example  6.3  was  given  by  More'  and  Cosnard  [10]  (also  see  More',  Garbow  and 
Hillstrom  [11]).  The  results  are  shown  in  the  tables  below,  where  IT  is  the  number  of  iterations, 
NF  is  the  number  of  function(F(x))  evaluations,  and  LN  is  the  number  of  line  searches  in  which 
the  step  length  X<1.  ND  is  the  number  of  nondecrease  directions.  xO  is  the  initial  guess. 


Example  6.1. 

/i(x)  =  8(xi-x|), 

Zy(x)  =  16x,(x/  -  Xy_,)  -  2(l  -  Xy)  +  8(xy-x/+1  ),  /=  2 .  »-l, 

/.(x)  =  16x,(x,2  -  x„_,)  -  2(1  -  x,), 
n  =  9  , 

xl  =  (-1,  -1 . -l)r,  x2  =  (-0.5,  -0.5,  ...,  -0.5)r,  x3  =  (2,  2,  ...,  2)r. 


Algorithms 

x0= 

=xl 

x0= 

=x2 

x0=x3 

IT 

NF 

LN 

ND 

IT 

NF 

LN 

ND 

IT 

NF 

LN 

ND 

CPR 

22 

88 

0 

22 

88 

15 

0 

8 

32 

■9 

0 

Schubert 

38 

41 

mm 

53 

56 

47 

5 

36 

ifia 

5 

Alg.  2.1 

fail 

'  1 

56 

114 

46 

14 

28 

Si 

0 

Alg.  5.1 

24 

50 

14 

mm 

24 

50 

15 

14 

30 

n 

0 

Table  6.1. 
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Example  6.2  (Broyden  tridiagonal  function). 

/.-(*)  =  (3  -  2 *,-)*,■  -  *,_!  -  2xf+1  +  1, 

zo  =  2»+i  =  0, 

n  =  9, 

*1  =  (-1,  -1,  -l)r,  x2  =  (-0.3,  0.3,  ...,  -0.3,  0.3)r, 

x3  =  (-10,  -10,  ...,  -10)r. 


Algorithms 

a:0= 

=xl 

x0==x2 

x0= 

=x3 

IT 

NF 

LN 

ND 

IT 

NF 

LN 

ND 

IT 

NF 

LN 

ND 

CPR 

5 

KM 

0 

0 

6 

24 

1 

0 

32 

0 

0 

Schubert 

7 

Ba 

0 

0 

11 

14 

2 

0 

3 

2 

Alg.  2.1 

8 

14 

0 

0 

8 

18 

2 

■ 

26 

0 

0 

Alg.  5.1 

6 

14 

0 

0 

7 

16 

2 

o 

11 

24 

0 

0 

Table  6.2. 


Example  6.S  (Discrete  boundary  value  function). 

l2 

fi(x)  =  2x,  -  z.,!  -  xf+1  +  — (  x,-  +  1,-  +  l)3 

A  =  1,-  =  ih,  x0  =  x,+1  =  0. 

n  =  9, 

zl  =  (»?/)r  ,  1j  —  -  1),  z2  =  (-1,  -1,  ...,  -l)r  , 

x3  =  (10,  10,  ...,  10)r. 


Algorithms 

x0= 

=xl 

x0=x2 

x0= 

=x3 

IT 

NF 

LN 

ND 

IT 

NF 

LN 

ND 

IT 

NF 

LN 

ND 

CPR 

3 

12 

0 

0 

4 

16 

0 

8 

32 

0 

0 

Schubert 

4 

mm 

5 

8 

0 

0 

17 

20 

2 

2 

Alg.  2.1 

4 

29 

0 

0 

6 

14 

12 

26 

0 

0 

Alft.  5.1 

4 

Bn 

0 

5 

12 

0 

o 

10 

22 

0 

0 

Table  6.3. 
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7.  Concluding  Remarks. 


We  have  presented  two  algorithms  for  solving  sparse  nonlinear  systems  of  equations.  The 
CM-successive  column  correction  algorithm  (Algorithm  2.1)  is  based  on  Coleman  and  More’s  par¬ 
titioning  algorithm  and  the  column-update  algorithm.  This  algorithm  uses  only  two  function 
values  at  each  iterative  step,  and  it  is  g-superlinearly  convergent.  Using  this  algorithm,  one 
group  of  the  columns  of  is  displaced  at  each  step.  Actually,  it  is  not  necessary  to  update  just 
one  group  at  each  iterative  step.  Instead,  we  can  displace  several  groups  at  each  iteration,  and 
this  gives  the  algorithm  a  faster  convergence  rate.  However,  if  one  more  group  is  displaced,  then 
one  more  function  value  is  required.  Therefore,  the  efficiency  of  the  algorithm  depends  on  the 
number  of  the  groups  displaced  at  each  iterative  step. 

The  modified  CM-successive  column  correction  algorithm  (Algorithm  5.1)  is  a  combination 
of  the  CM-successive  column  correction  algorithm  and  Schubert’s  algorithm.  It  is  also  q- 
superlinearly  convergent.  Our  numerical  results  indicate  that  the  modified  successive  column 
correction  algorithm  usually  behaves  much  better  than  the  CM-successive  column  correction  algo¬ 
rithm.  However,  we  have  not  been  able  to  prove  better  theoretical  convergence  results  for  the 
modified  CM-successive  column  correction  algorithm  than  those  for  the  unmodified  one.  Addi¬ 
tional  numerical  results  indicate  that  the  modified  CM-successive  column  correction  algorithm  is 
also  usually  more  efficient  than  the  CPR-CM  algorithm  and  Schubert’s  algorithm.  When  the  prob¬ 
lem  is  not  well  behaved,  or  the  initial  guess  is  far  away  from  the  solution,  the  modified  CM- 
successive  column  correction  algorithm  is  much  more  efficient  than  Schubert’s  algorithm. 

The  idea  of  the  CM-successive  column  correction  algorithms  can  also  be  used  with  Powell 
and  Toint’s  [14]  work,  which  will  lead  to  methods  for  unconstrained  optimization  problems.  This 
will  be  our  future  work. 
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